COMPLETE SUBDIVISION ALGORITHMS, II: 
ISOTOPIC MESHING OF SINGULAR ALGEBRAIC CURVES 



MICHAEL BURR, SUNG WOO CHOI, BEN GALEHOUSE, AND CHEE YAP 



Abstract. Given a real valued function f(X,Y), a box region S C K 2 and 
e > 0, we want to compute an e-isotopic polygonal approximation to the re- 
striction of the curve S = / _1 (0) = {p £ K 2 : f(p) = 0} to B . We focus on 
subdivision algorithms because of their adaptive complexity and ease of imple- 
mentation. Plantinga & Vegter gave a numerical subdivision algorithm that 
is exact when the curve S is bounded and non-singular. They used a compu- 
tational model that relied only on function evaluation and interval arithmetic. 
We generalize their algorithm to any bounded (but possibly non-simply con- 
nected) region that docs not contain singularities of S. With this generalization 
as a subroutine, we provide a method to detect isolated algebraic singularities 
and their branching degree. This appears to be the first complete purely nu- 
merical method to compute isotopic approximations of algebraic curves with 
isolated singularities. 

Key words: Meshing, Singularity, Root bound, Evaluation bound, Implicit 
algebraic curve, Complete numerical algorithm, Subdivision algorithm. 



1. Introduction 

Given s > 0, a box region Bq C ]R 2 and a real valued function / : W 2 — > K, 
we want to compute a polygonal approximation P to the restriction of the implicit 
curve S = / _1 (0) to B (where / -1 (0) = {p £ R 2 : f(p) = 0}). The approximation 
P must be (1) "topologically correct" and (2) "e-close" to S fl B . We use the 
standard interpretation of requirement (2), that d(P,SP\Bo) < e where d(-,-) is 
the Hausdorff distance on compact sets. In recent years, it has become accepted 



Boissonnat et al. 2006 to interpret requirement (1) to mean that P is isotopic to 
S n Bq , which we denote by P ~ S fl Bo ■ This means that we not only require 
that P and S n B are homeomorphic, but also require that they are embedded in 
M 2 "in the same way" . This means that the two embeddings can be continuously 
deformed to each other. E.g., if S Pi Bo consists of two disjoint ovals, these can 
be embedded in M 2 as two ovals exterior to each other, or as two nested ovals. 
Isotopy, but not homeomorphism, requires P to respect this distinction. There is 
a stronger notion of isotopy called ambient isotopy (see the definition in Section 
4). We use this stronger notion in this paper (but, for simplicity, we still say 



"isotopy"). See Boissonnat et al. 2006 p. 183] for a discussion of the connection 



between ambient and plain isotopy. In this paper, we focus mainly on topological 
correctness since achieving £-closeness is not an issue for our particular subdivision 
approach (cf. |Boissonnat et al." 2006 pp. 213-4]). This amounts to setting e = oo. 
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We may call the preceding problem the 2-D implicit meshing problem. The 

term "meshing" comes from the corresponding problem in 3-D: given e > and an 
implicit surface S : f(X, Y, Z) = 0, we want to construct a triangular mesh M such 
that d(M, S) < e and M ss S. It is interesting (see Burr et a3T| |2009b] ) to identify 



the 1-D meshing with the well-known problem of real root isolation and refinement 
for a real function f(X). 

The algebraic approach and the numerical approach constitute two ex- 
tremes of a spectrum among the approaches to most computational problems on 
curves and surfaces. Algebraic methods can clearly solve most problems in this 
area, e.g., by an application of the general theory of cylindrical algebraic decompo- 
sition (CAD) |Basu et aL 2003 . Purely algebraic methods, however, are generally 



not considered practical, even in the plane (e.g., Hong 1996 Seidel and Wolpert 
2005] ), but efficient solutions have been achieved for special cases such as inter- 



secting quadrics in 3-D Schoemer and Wolpert 2006 . At the other end of the 
spectrum, the numerical approaches emphasize approximation and iteration. An 
important class of such algorithms is the class of subdivision algorithms which 
can be viewed as a generalization of binary search. Such algorithms are practical in 
two senses: they are easy to implement and their complexity is more adaptive with 
respect to the input instance |Yap| |2006| . Another key feature of subdivision algo- 
rithms is that they are "localized", meaning that we can restrict our computation 
to some region of interest. 

Besides the algebraic and numerical approaches, there is another approach that 
might be called the geometric approach in which we postulate an abstract com- 
putational model with certain (geometric) primitives (e.g., shoot an ray or decide 
if a point is in a cell). When implementing these geometric algorithms, one must 
still choose an algebraic or numeric implementation of these primitives. Implemen- 
tations can also use a hybrid of algebraic and numeric techniques. 

Unfortunately, numerical methods seldom have global correctness guarantees. 



The most famous example is the Marching Cube algorithm Lorensen and Cline 



1987 



Many authors have tried to improve the correctness of subdivision algorithms 
1997] ). So far, such efforts have succeeded under one of 



Standcr and Hart 



(e.g., 

the following situations: 

• (AO) Requiring niceness assumptions such as being non-singular or Morse. 

• (Al) Invoking algebraic techniques such as resultant computations or ma- 
nipulations of algebraic numbers. 

It is clear that (AO) should be avoided. Generally, we call a method "complete" 
if the method is correct without any (AO) type restrictions. But many incomplete 
algorithms (e.g., Marching cube) are quite useful in practice. We want to avoid 
(Al) conditions because algebraic manipulations are harder to implement and such 
techniques are relatively expensive and non-adaptive 



Yap] 12006] . The complete 



removal of (AO) type restrictions is the major open problem faced by purely nu- 
merical approaches to meshing. Thus, Boissonnat et al. 2006 p. 187] state that 



"meshing in the vicinity of singularities is a difficult open problem and an active 
area of research" . Most of the techniques described in their survey are unable to 
handle singularities. It should be evident that this open problem has an implicit 
requirement to avoid the use of (Al) techniques. 



For example, the subdivision meshing algorithm of Plantinga 2006 Plantinga 



and Vegter 2004 requires the non-singularity and the boundedness of curves and 
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surfaces, (AO) assumptions. The subdivision algorithm of Seidel and Wolpert 
require the computation of resultants, an (Al) technique. We thus classify 



2005 



Seidel 



and Wolpert| 2005 as a hybrid approach that combines numerical and algebraic 
techniques. Prior to our work, we are not aware of any meshing algorithm that 
can handle singularities without resorting to resultant computations. In general, 
hybrid methods offer considerable promise (e.g., [Hong [19 96 ). This is part of a 
growing trend to employ numerical techniques to speed up algebraic computations. 



Some of our recent work addresses the above (AO) / (Al) concerns: in Yap 2006 



we gave a complete numerical approach for determining tangential Bezier curve in- 

we numerically solve zero-dimensional triangu- 



Chcng et al. 2007 



tersections; in 

lar systems without any "regularity" requirements on the systems; in Burr et al. 



2009b , we provide numerical root isolation in the presence of multiple zeros; and 



Burr et al. 2009a provides one of the first non-probabilistic adaptive analyses of 



an evaluation-based real root isolation algorithm. These last two papers address 
the 1-D analogue of the Plantinga & Vegter Algorithm. The philosophy behind all 
of these papers is the design and analysis of complete numerical methods based on 
approximations, iteration and adaptive methods. Topological exactness is achieved 
using suitable algebraic bounds, ranging from classical root separation bounds to 
evaluation bounds and geometric separation bounds. We stress that the worst-case 
complexity of adaptive algorithms ought not to be the chief criterion for evaluat- 
ing the usefulness of these algorithms: for the majority of inputs, these algorithms 
terminate fast. Zero bounds are only used as stopping criteria for iteration in the 
algorithms, and simple estimates for them can be computed easily. Computing such 
bounds does not mean that we compute resultants, even though their justification 
depends on resultant theory. The present paper continues this line of investigation. 
The recent collection Boissonnat et al. 2006, Chapter 5] reviews the current 



algorithmic literature in meshing in 2- and 3-D: the subdivision approach is rep- 
resented by the Plantinga & Vegter algorithm as well as by Snyder's earlier ap- 
proach based on parametrizability Snyder 1992b|a . The subdivision algorithm of 
Plantinga & Vegter is remarkable in the following sense: even though it is glob- 
ally isotopic, it does not guarantee isotopy of the curve within each cell of the 
subdivision. In contrast, Snyder's subdivision approach Snyder 1992b|a requires 
the correct isotopy type in each cell. Indeed, because of this, the algorithm is 
incomplete [Boissonnat et al. 1 120061 p. 195]. 



Among geometric approaches to meshing, we have the point sampling approach 



as represented by Boissonnat and Oudot 2006 Cheng et al. 2004 , the Morse The- 



ory approach as represented by Stander and Hart 1997 Boissonnat et al. 2004 



and the sweepline approach Mourrain and Tecourtl 2005 . Note that the sweepline 



approach naturally corresponds to the algebraic operation of projection; therefore, 
its implementation is often purely algebraic. The idea of the sampling approach 
is to reduce meshing of a surface S to computing the Delaunay triangulation of a 
sufficiently dense set of sample points on S Boissonnat et al. 2006 p. 201-213]. 



To obtain such sample points, Cheng et al. 2004 need a primitive operation that 



amounts to solving a system of equations involving / and its derivatives. Boisson- 



nat and Oudot 2006 need a primitive for intersecting the surface with a Voronoi 



edge. These sample points are algebraic, and implementing the primitives exactly 



1 Their paper is subtitled "Exploiting a little more Geometry and a little less Algebra" which 
speaks to our concerns with (Al). 
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would require strong algebraic techniques. But exact implementation docs not seem 
to be justified for these applications, and so we are faced with an implementation 
gap that shows well-known non-robustness issues. For restrictions and open prob- 



lems in sampling approaches, see Boissonnat et al. 2006 p. 227-229]. In contrast 



the computational primitives needed by subdivision approachs work directly with 
bigfloats, with modest requirements on /. 

This paper presents a purely numerical subdivision method for meshing alge- 
braic curves with isolated singularities. In a certain sense, this is the most general 
geometric situation because, by Proposition [l] reduced algebraic curves have only 
isolated singularities. Our starting point is the algorithm of Plantinga & Vegter 
Plantinga and Vegter, 2004 Plantinga 2006 for implicit meshing of curves. It is 



important to understand the computational model of Plantinga & Vegter which is 
also used in this paper. Two capabilities are assumed with regards to f(X,Y): 

• (i) Sign evaluation of f(p) at dyadic points p. 

• (ii) / is C 1 and we can evaluate the interval analogues (i.e., box functions) 
oi f, §y on dyadic intervals. 

Note that the Marching Cube algorithm only requires capability (i). Let the class 
PV denote the set of all real functions / : M 2 — > R for which capabilities (i) and (ii) 
are available. Many common functions of analysis (e.g., the elementary functions 



Du and Yap 2006| ) belong to PV; thus, the approach of Plantinga & Vegter admits 



a more general setting than algebraic curves. 

1.1. Overview of Paper. 

• Section 2 establishes some basic terminology and recalls facts about the 
singularities of algebraic sets. 

• In Sections 3 and 4, we extend the Plantinga & Vegter algorithm to compute 
an isotopic approximation of the curve S = / _1 (0) restricted to a "nice 
region" that need not be simply connected. S may have singularities outside 
R and we only need f £ PV. 

• In Section 5, we provide the algebraic evaluation bounds necessary for mesh- 
ing singular curves. 

• In Section 6, we provide a subdivision method to isolate all the singularities 
of a square-free integer polynomial f(X,Y), 

• In Section 7, given a box B containing an isolated singularity p, we provide 
a method to compute the branching degree of p. 

• In Section 8, we finally present the overall algorithm to compute the isotopic 
polygonal approximation. 

• We conclude in Section 9. 



2. Basic Terminology and Algebraic Facts 

Let F:=Z [f] = {m2™ : m,n € Z} be the set of dyadic numbers. All of our 
numerical computations are performed in a straightforward manner using F. There 
are many well-known implementations of arithmetic on such numbers, and, in this 
case, they known as bigfloats. In short, our computational model is not based 
on some abstract capability whose implementation may reveal gaps that lead to 
well-known non-robustness issues. 

For S C R, let US be the set of closed intervals [a, b] with endpoints in S,a,b € S. 
We write US n for (US) n . In particular, fJF is the set of dyadic intervals, and UM. n is 
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the set of n-boxes. The width of an interval I — [a, b] is w(I) :=b — a. The width 
and diameter, respectively, of an n-box B — n"=i[ a *'^] * s W (B) :=min{6j — : 
i = 1, . . . ,n} and d(B) :=max{&i — : i — 1, . . . ,n}. The boundary of a set 
SCRis denoted dS. If / : K™ — >• K and SCI, then f(S) :={f(x) : x G S}. 
A function □/ : OF™ ->■ DF is a box function for / provided (i) f (B) C Df(B) 
and (ii) if Bo D B\ D ■■ ■ with lim^ Bi — p then linijD = f(p). We regard 
the limit of intervals in terms of the limits of their endpoints. We say / € PV 
if / € C 1 (has continuous first derivatives), there is an algorithm to determine 
sign(/(p)) for p € F™ and □/ and the corresponding functions for the derivatives 
of / arc computable in DF. In this paper, we only consider box functions for the 
two dimensional case. 

We only consider boxes of the form B = I x J where /, J are dyadic intervals. 
As in Plantinga & Vegter, all our boxes will be squares, w(I) = w(J) = w(B). Our 
algorithms work in the slightly more general setting where all boxes have aspect 



ratios at most 2 (see Lin and Yap 2009 for a proof). The boundary dB of B is 



divided into four sides and four corners. Note that the 'sides/corners' terminology 
for boxes should not be confused with the 'edges/vertices' terminology which we 
reserve for the straightlinc graph G = (V, E) which is the approximation to our 
curve. We split a box B by subdividing it into 4 subboxes of equal widths. These 
subboxes are the children of B and each has width ^w(B). Starting with Bo, 
the child-parent relationships obtained by an arbitrary sequence of splits yields a 
quadtree rooted at Bo- Two distinct boxes B, B' of a quadtree are neighbors if 
their boundaries overlap: B C\B' is a line segment, but not a single point. Segments 
of the form BOB' are called interior segments. If B is a box whose side s is part 
of the boundary of Bo, then we call s a boundary segment and B a boundary 
box. Moreover, each side of a box is divided into one or more segments. 

Some of the regions of interest in this paper will not be squares and may not even 
be simply connected. To ensure that our algorithm continues to work in this case, 
we insist that the region comes from a subdivision. I.e., we insist that there exists 
a square Bo, a subdivision of Bo and a collection of boxes in the subdivision such 
that the region of interest is the union of this collection of boxes. Although it is 
not necessary for our algorithms, in most implementations, it is easiest to maintain 
a subdivision of Bq where each box is labeled as either region or complement. The 
union of the boxes labeled "region" will be written R (or Rq(T) when we stress 
the subdivision tree) and is the region of interest. We reserve the notation Ro for 
non-square regions and use Bo for square regions. 

Algebraic Facts. Let D be a unique factorization domain (UFD) and f,g£ 
D[X] = D[Ai, . . . , X n ] where X = (Xi, . . . , X n ). We say /, g are similar if there 
exist a, b <E D \ {0} such that af = bg, and we write this relationship as / ~ g. 
Otherwise, / and g are dissimilar. The square-free part of / is defined as 

(1) SqFree(/) - 



GCD(f,d Xl f,...,d x J) 



where dxi indicates differentiation with respect to Xi. f is said to be square- free 
if SqFree(/) = /. From |l]) we see that computing SqFree(J) from / involves 
only rational operations of D. As the gradient of / is V/ = (c?Xi/, • ■ • , dx„f), we 
may also write GCD(/, V/) for GCD(/, d Xl f,..., d x J). See |Yap][2llu0l Chap. 2] for 
standard conventions concerning GCD. 
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Let k be an algebraically closed field. For S C fc[X] = k[X\, . . . ,X n ], let 
Zero(S') :—{p £ k n : f(p) = for all / £ S} denote the zero set of S. A zero set 
is also known as a variety. The singular points of Zero(/) are defined to be the 
points where VSqFree(J) = 0. 

In 1-dimension, it is well-known that a square-free polynomial / G Z[X] has no 
singularities (i.e., multiple zeros). We now recall two generalizations of this result 
that will be necessary in the remainder of the paper. See |Hartshorne| |1977| Cox 



et al. . 1992 Harris 1992 for similar results 



PROPOSITION 1 ( [Harris 1992 Ex. 14. 3]). The singular points of any variety form 



a proper subvariety. 

This proposition is critical in our paper, because it implies that if / £ M[X, Y] 
is square-free, then the singular points are a proper subvariety of a union of curves 
and hence must be a finite set of points. Thus, we only need to handle isolated 
singularities. 



Proposition 2 (Algebraic Sard Lemma Harris 1992 Prop. 14.4]). Let f : X -> Y 



be any surjective regular map of varieties defined over a field k of characteristic 0. 
Then there exists a nonempty open subset U C Y such that for any smooth point 
p £ l~l X sm in the inverse image ofU, the differential df p is surjective. 

Here, X sm denotes the set of smooth points of variety X. The open sets refer 
to the Zariski topology. The condition that the differential df p is surjective is 
equivalent to insisting that the Jacobian of / has the same rank as the dimension 
of Y. The situation that we consider is / : R 2 — ► R where / is a polynomial. In this 
case, since the image is one dimensional, the condition that df p is surjective reduces 
to the condition that V/(p) 7^ 0. Every point in R 2 = X is smooth and R \ U is 
only a finite set. Hence, there are only a finite number of level sets, parametrized 
by h, where Zero(/(A, Y) — h) has a singular point. 

3. Algorithm of Plantinga & Vegter 

First, we recall the Plantinga & Vegter algorithm: Given s > 0, a bounded and 
nonsingular curve / _1 (0) for / : R 2 — >• R and a bounding box Bo £ OF 2 , they 
compute a polygonal curve P. The correctness statement is as follows: if f~ 1 (0) Q 
Bq, then P is an e- approximation to the curve S = / (0), i.e., d(P,S) < e and 
P fa S. For simplicity, they focus on topological correctness: P ps S, since it is easy 
to refine the subdivision to achieve d(P, S) < e. The curves in Figure [l] are output 
from our implementation of the Plantinga & Vegter algorithm as reported i n |Lin| 
and Yap , 2009 . The value of e is small in Figure [IJa) , while e — 00 in Figure |ljb) . 



The algorithm is based on two simple predicates on boxes B: 

• Predicate C {B) holds if ^Uf(B). 

• Predicate d(B) holds if g (□ J£(B)) + 

These predicates are easily implemented for / £ PV. If Cq(B) holds, then the 
curve S does not intersect B. If C\(B) holds, then the gradient of / is never zero in 
B and the gradient vectors approximately point in the same direction. Note that if 
B satisfies Ci, then, by recalling the parent-child relationship, any child of B also 
satisfies C\. 
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The input box B is a dyadic square and the output will be an undirected graph 
G = (V, E) where each vertex v <E V is a dyadic point, v G F 2 . In fact, G represents 
a straightline planar graph that is a polygonal approximation of S. 

The algorithm has 3 phases, where Phase i [i = 1, 2, 3) is associated with a queue 
Qi containing boxes. Initially, Q\ = {-Bo}, an d Qi = Qz = 0- When Qi is empty, 
proceed to Phase i + 1. 

• PHASE 1: SUBDIVISION. While Qi is non-empty, remove some B from 
Qi, and perform the following: If Cq(B) holds, discard B. If C\(B) holds, 
insert B into Qi. Otherwise, split B into four subboxes and insert them 
into Qi. 

• PHASE 2: BALANCING. This phase "balances" the subdivision, where a 
subdivision is balanced if the widths of any two neighboring boxes differ 
by at most a factor of 2. Queue Qi is a min-priority queue, where the 
width of a box serves as its priority. While Qi is non-empty, remove the 
min-element B from Qi, and perform the following: For each S-neighbor 
B' in Qi with width more than twice the width of B, remove B' from Q 2 
and split B' . Insert each child B" of B' into Qi provided C (B") does not 
hold. B" might be a new neighbor of B and B" might be split subsequently. 
Finally, when every neighbor of B is at most twice the width of B, insert 
B into Q3. 

• PHASE 3: CONSTRUCTION. This phase constructs the graph G = (V, E). 
Initially, the boxes in Q 3 are unmarked. While Q3 is non-empty, remove 
any B from Q3 and mark it. Now construct the set V(B) of its vertices. 
For each B-neighbor B' , if B' is marked, retrieve any vertex v on the side 
B n B' , and put v into V(B). If B' is unmarked, evaluate the sign of 
f{p)f(q) where p, q are endpoints of the segment B n B' . If f(p)f(q) < 0, 
create a vertex v — (p + q)/2 for the graph G and put v into V(B). Note 
that if f(p) = for any corner p, treat f(p) as positive; in effect, this is an 
inhnitesimal perturbation at p. It can be shown that |V(i?)| €E {0,2,4}. If 
|V(I?)| = 2, put the edge (p, q) into G to connect the vertices in V(B). If 
|U(i?)| = 4, it can be shown that one side of B contains two vertices a, b. 
Introduce two edges into E to connect each of a, b to the remaining two 
vertices. The requirement that these two edges are non-crossing ensures 
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that the connection is unique (see Plantinga and Vegter 2004 Plantinga 
20061 ). 



The output graph G — (V, E) can be viewed as a straightline graph, decomposed 
into a collection P = P(G) of closed polygons. In the following, we simply use G 
in place of P as the polygonal approximation. 

4. Extension of Plantinga & Vegter 

4.1. The Naive Extension of Plantinga & Vegter. The correctness statement 
of Plantinga & Vegter requires Bq to be a bounding box for the curve / _1 (0). The 
power of subdivision methods comes from their ability to adaptively analyze local 
data. Choosing the initial box to be a bounding box gives up this power of local 
analysis; also, by insisting that the curve is bounded, the types of possible curves 
are severely restricted. As a first attempt, one could naively attempt to run the 
Plantinga & Vegter Algorithm starting with an arbitrary box and then ask the 
question "In what sense is the output G correct?" Intuitively, G should be isotopic 
to / _1 (0) l~l Bo, but Plantinga & Vegter did not discuss this issue. The algorithm 
certainly cannot handle the case when the curve S = f~ 1 (0) has tangential but 



non-crossing intersections Yap 2006 with 8Bq. If we assume that there are only 
transversal intersections, we still face two problems: if the curve S (locally) enters 
and exits OBq by visiting only one box B C Bo, the above algorithm would fail to 
detect this small component. Such an error is called an undetected incursion, as 
illustrated in Figure J2^a). Conversely, the curve S might escape undetected from 
Bq. Such an error is called an undetected excursion, as illustrated in Figure [2jb). 
These errors cause the final approximation to be incorrect, since S P\ B may not 




Figure 2. (a) incursion, (b) excursion, (c) boundary boxes and 
their complements 



have the same number of components as G. Incursions account for components of 
S n B that are undetected by the Plantinga & Vegter algorithm, and excursions 
account for connected components of G that are not connected in SdB . A single 
box may contain several incursions and excursions. The excursion case is more 
troubling because there is no guarantee that C\ will hold in the complement of 
Bq; this means that the connected components of G which are not connected in 
S n B might, in fact, not even be connected in S. If we choose B large enough, 
such errors cannot arise, but this approach gives up the power of adaptivity and 
localization which subdivision methods possess. If S has singularities, making Bq 
large may not be an option. In this paper, we avoid any "largeness" assumption on 

The preceding issues arise because the Plantinga & Vegter algorithm focuses 
only on the parity of the endpoints of an edge of the subdivision. They prove 
that multiple intersections can be removed by applying a suitable isotopy to move 
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small features into a neighboring box. We reproduce their results here in order to 
examine them. In the following, we fix a balanced quadtree T rooted at a box Bq. 
Let V(T) denote the collection of boxes at the leaves of T. These boxes constitute 
a partition of Bq. We also assume that Cq(B) or C\(B) holds at each B £ V(T). 
The subdivision constructed through the Plantinga & Vegter algorithm has this 
property. The results in this section do not need the tree T to be balanced because 
this constraint is only used in the construction of the isotopic curve. 

In PHASE 3 of the Plantinga & Vegter algorithm, we noted that if / _1 (0) passes 
though a corner of the subdivision, then they treat that point as positive. This 
can be done via an infinitesimal perturbation of /: we call this / the standard 
perturbation of / with respect to T. More precisely, / is defined to be a function 
that agrees with / everywhere except in infinitesimal neighborhoods of corners of 
boxes where / is 0; in such neighborhoods, / is positive at the corner of the box. 
Clearly, the definition of / depends on T. In general, S = / _1 (0) is isotopic to 
S = / _1 (0), but in the case where Bq is not a bounding box for the curve, it can 
happen that S (1 Bp is n ot isotopic to S PI Bq. This difference is discussed more 
throughly in Section 



4.3 



All of our correctness statements are about S and not S. 
However, for the sake of simplicity in this section, we continue to use the notations 
S and / instead of S and /. 

Consider the dual graph H = H(T) whose vertex set is V(T) and edges connect 
pairs of neighboring boxes. Thus, B, B' € V(T) are connected by an edge in H 
iff s = B PI B' is a segment, but not a point. We are interested in the subgraph 
Hf = Hf(T) of H in which B,B' are connected by an edge iff the segment s — 
BOB' intersects the curve / _1 (0) more than once. We further focus on the non- 
trivial connected components of Hf, i.e., those components consisting of more than 
one box. In Figure [3j the graph Hf has three non-trivial connected components 
(A,B,C). The union of all of the boxes corresponding to a non-trivial connected 
component of Hf is called a monotone region, and Hf itself will be known as the 
monotone graph. In Figure |3][a) , we have filled each monotone region (A,B,C) 
with a distinct tint. 




FIGURE 3. A quadtree T and its dual graph induced by / 1 (0). 



Let us explain the monotone terminology. Each monotone region M can be 
uniquely classified as a V-monotone or A-monotone region in the following sense: 
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we say M is F-monotone relative to / if for all B,B' C M where (B,B') is 
an edge of Hf, then B n B' is a segment orthogonal to the ^-direction (i.e., s is 
horizontal). Moreover, 

(2) 0£^(BUB>). 

We call such a rectangular region BUB' a local neighborhood of the monotone 
region. Thus ^ implies that the curve / _1 (0) restricted to each local neighborhood 
of a y-monotone region is parametrizablc in the A-direction (i.e., each vertical line 
intersects the curve at most once). This is a "local property" because the curve 
/ _1 (0) might not be parametrizable in the X-direction in the entire K-monotone 
region. See Figure [4] for an example of a X-monotone region in which the curve 
/ _1 (0) is not globally parametrizable in the Indirection: a horizontal line L inter- 
sets the curve twice in the monotone region. The curve is parametrizable, however, 
for connected components of L with respect to the X-monotone region. The defi- 
nition of a X-monotone regions is similar, after exchanging the roles of X & Y and 
the roles of horizontal & vertical. 




Figure 4. An X-monotone region that is not globally parametriz- 
able in the ^-direction. 



Let C be a connected component of / _1 (0) n B. We classify C into three types: 
if both endpoints of C lie on one side of B, then C is an incursion as discussed 
above. If the endpoints of C lie on two adjacent sides of B, it is called a corner 
component. Finally, if the endpoints of C lie on two opposite sides of B, it is 
called a crossing component. These three kinds of components are illustrated 
in Figure [5^ a-c) . Note that there can be no other types of components in B. In 
particular, B cannot contain an isolated component since C\{B) holds. The next 
lemma also restricts the number of crossing components. 

Let us define the notion of isotopy used in this paper. An isotopy is a continuous 
function I : R 2 x [0, 1] — » K 2 such that the following two conditions hold: (1) /(•, 0) 
is the identity and (2) /(•, t) is a homeomorphism for each t € [0, 1]. If, in addition, 
there is a set BCI 2 such that for all t £ [0, 1], I(-,t) acts as the identity on the 
complement of B, i.e., I(p, t) = p for allp ^ B, then we call / an isotopy on B, and 
write the isotopy as I : B x [0, 1] —> B. Two sets S, S' C K are said to be isotopic, 
denoted S w S' , if there exists an isotopy I such that I(S, 1) = S'. As promised 



in the introduction, our notion of isotopy is the "ambient" sort Boissonnat et al. 



2006 p. 183]. Observe that if I is an isotopy on B, then for any continuous function 
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Figure 5. Types of components and their isotopic transformation 
in a local neighborhood. 



/:£?—>• R, / _1 (0) and (/o/- 1 ^,*))" 1 ^) are isotopic. We state without proof the 



following lemma which is somewhat implicit Plantinga and Vegter 2004 Plantinga 



2006 



Lemma 3. Let (B,B') be an edge of the monotone graph Hf. Assume their shared 
segment s = B l~l B' is horizontal. 

(1) B U B' is Y -monotone relative to f 

(2) There is at most one crossing component in each of B and B' . If there is 
one crossing component in each of B and B' then these two components 
are part of one single crossing component of B U B' . 

(3) Lef/^ 1 (0) intersect s at two consecutive points p\ andp2- Then pi andpi 
are connected by a curve component X{p\,p2) C / _1 (0) that lies entirely 
within B or entirely within B' . 

(4) There is an isotopy I on B U B' that reduces the number of intersections 
(counted with multiplicities) on s so that (/ o /(•, l)" 1 )" 1 ^) intersects s 
in 2 fewer times than / _1 (0), and B U B' remains Y -monotone relative to 
/o/(-, t )- 1 for allt. 

Part (2) of this Lemma follows from the fact that Ci(BUB') holds. The situation 
where B and B' both has a crossing component in illustrated in Figure [5^(1) . The 
isotopy lonBUB' constructed in part (4) cannot decrease the number of crossing 
components, but may increase the number by one, as illustrated in Figure^e). The 
general idea is to repeatedly apply the transformation given by such an isotopy / 
until the curve / _1 (0) composed with the inverses of the isotopies intersects the 

for a generalization under a weaker 



2009 



segment s at most once (see Lin and Yap 
predicate than C\). The isotopies can easily be chosen to preserve monotonicity 
by keeping the appropriate coordinate fixed throughout the isotopy, i.e., the X- 
coordinate can be kept fixed for an isotopy on a y-monotone regions. This prepares 
us for an induction because after each isotopic transformation, monotonicity is 
preserved. To develop the appropriate global correctness statement, we formulate 
the sense in which these isotopies interact: 

Lemma 4. Let H be the dual graph to the subdivision given by the tree T , i.e., there 
is an edge connecting boxes B\ and B2 iff they are neighbors. Let Hf be the subgraph 
of H where B\ and Bi are connected by an edge of H and / _1 (0) intersects their 
shared segment more than once. Let D be a connected component of Hf, and R be 
the union of the boxes appearing in D . There exists an isotopy I on R such that if 
fi = f then f- 1 (Q)nR^ fr~ 1 (0)nR andf^iO) intersects all segments 

between boxes of D at most once. 
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Note, in particular, that boxes which correspond to isolated components of Hf 
are boxes where / _1 (0) intersects each of its sides at most once. The upshot of this 
argument may be summarized by the following theorem, implicit in Plantinga & 
Vegter: 

Theorem 5. Let T be the subdivision tree at the end of the Plantinga & Vegter 
algorithm. Then there exists an isotopy I on Bq such that / _1 (0) n -Bo ~ fx (0) H 
B , where f\ := / o /(•, Moreover, /^ _1 (0) intersects each interior segment of 

T at most once. 

A function f\ as given by the conditions of this theorem, is said to be normalized 
relative to T. It is important to realize that the normalized curve / 1 ~ 1 (0) may 
intersect the boundary segments of T more than once. The reason for this is that 
we can only apply the transformation of Lemma [3] to interior segments; we cannot 
apply the transformation to boundary segments without changing the topology 
of the restriction of the curve to Bq. Therefore, since Plantinga & Vegter focus 
entirely on the parity across segments, the graph G constructed by the Plantinga 
& Vegter algorithm may not be isotopic to /]" (0) n Bq. We note that although 
the isotopies constructed here might not preserve the C\ condition since they are 
only guaranteed to preserve monotonicity, it will be useful to recall that these boxes 
were derived from ones where C\ did hold. 

One solution for ensuring that G satisfies G f=s S fl Bq is to ensure that S 
intersects each boundary segment at most once. A sufficient condition is that 
either Cq or C\ holds on the edge, using the corresponding properties for the one- 
dimensional version of Plantinga & Vegter 's algorithm (cf. [Burr et al. 2009b] ). 



Achieving this sufficiency condition requires that we determine the topology of 
S on the boundary of Bq, including tangental intersections. This solution is also 
inefficient because it may require frequent subdivisions near the boundary to resolve 
fine features; in higher dimensions, this solution must also be recursively applied 



to lower dimensional boxes. These issues also arise in Snyder's approach Snyder 



1992b 



4.2. The Enlarged Region Solution. We now provide an alternative solution 
which is closer to the spirit of the Plantinga & Vegter approach of exploiting isotopy. 
This solution avoids determining the exact boundary topology as well as making 
any largness asssumptions on Bo. We wish to slightly enlarge Bq so that incursions 
and excursions can be removed by an appropriate isotopy. The basic idea is that, in 
addition to subdividing Bq , we find a slightly larger region B' which is the union of 
Bq with a "collar" of squares around Bq. We construct a straightline approximation 
G + that will be isotopic to the curve restricted to some expansion Bq of Bo into 
this collar. This is weaker than saying G + is isotopic to / _1 (0) fl Bq, but it has 
two major advantages: it allows our algorithm to terminate faster and it does not 
require / _1 (0) to intersect 8Bq in any special way. 

Call a box B C Bq a boundary box if dB intersects dB . Let B be such a box. 
If B does not share a corner with Bq, then it has a unique complementary box 
B' such that B' has the same width as B, the interiors of B' and Bq are disjoint, 
and dB' n 8Bq = dB n 8Bq. B and B' are called partners of each other. If B 
shares a corner with Bq, then it determines two complementary boxes B' ,B" . See 
Figure [2jc) . We insist that the collar is formed from complementary boxes all of 
which satisfy either C or C\. It is easy to adapt Plantinga & Vegter's algorithm 
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so that this situation occurs. With this property, excursions from Bq or incursions 
into B can be limited to to the collar region. 

Let B be a boundary box with complementary box B' . If the complementary 
box satisfies Co, then the curve does not intersect it. Therefore, there can be no 
incursions into or excursions from B. If the complementary box satisfies C\ and 
there is an excursion, then the isotopy presented in Plantinga & Vegter shows that 
the approximation constructed by the naive Plantinga & Vegter algorithm is correct 
in B U B l+ , for some region B C B' + C B' U B. The case that is left to consider is 
when there is an incursion. This case is harder because the isotopy constructed by 
Plantinga & Vegter would remove the component of the incursion, but that would 
result in an error. In the spirit of Plantinga & Vegter's algorithm, we consider the 
sign pattern of the corners of complementary boxes. 




Figure 6. Classification of complementary boxes according to 
sign of / at the corners. 



Up to symmetry and sign flips, B' can be put into one of five types based solely 
on the sign of / at the corners of B' , (a)-(e). This is illustrated in Figure [6j Note 
that the "alternate sign pattern" does not appear because this pattern cannot be 
C\ as shown by Plantinga & Vegter. 

We draw two instances of Types (a) to indicate the two possible "dispositions" 
of the curve / _1 (0) in B': In (al) the curve makes no incursion into B, but (a2) 
represents at least one incursion into B. These two dispositions give rise to dis- 
tinct isotopy types for the curve / _1 (0) n B. Similarly, type (b) has two possible 
dispositions (but we only indicate the case where there is an incursion). Types 
(c), (d), and (e) do not have analogous dispositions. Because of this difference in 
dispositions, we further classify Types (a) and (b) as transient; the other types 
are called non-transient. 

Some Intuition. We show how the complementary boxes are used to yield a 
correctness statement. Suppose B' is a complementary box whose partner is B. 
Let the straightline G = (V,E), when restricted to B, be denoted GC\B. Note that 
GOB has at most two edges. We would like to claim that /~ 1 (0) n B is isotopic to 
GOB. This is evidently false for Types (a) and (b) because of the two dispositions 
discussed above; but even for Type (c), this claim can be false because / _1 (0) n B 
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may have an arbitrary number of components due to incursions/excursions, as 
illustrated in Figure W[c) . 



B 



o 



B B 

— t r — — ? — — ? 

(a) Extending B by a "cap" B' + (b) / _1 (0) fl B + (c) Isotopic approximation 




Figure 7. Representing the transient type: extension of a bound- 
ary box by a cap 



To remedy the situation for Types (a) and (b), we first expand B into a slightly 
larger region B + C B U B' . The actual expansion B' + := B + \ B is represented 
in Figures [6| ^\&,h) by the white region demarcated by a dashed curve and we 
call B' + a cap. Up to isotopy, the incursions into B + in Types (al) and (a2) are 
equivalent to a single incursion of the curve / _1 (0) into B + from B' , as shown in 
Figure [7|h). In the graph G — (V, £7), we represent this incursion by adding three 
vertices u, v, w to y, with u and u> on the side B n £? , and with v between m and 
w, but slightly to the interior of B. We also add the edges (u, v), (v, w) to E. This 
is illustrated in Figure [T^c). Let G + denote the augmented straightline graph 
with these additional vertices and edges for transient complementary boxes. 

Types (a) and (6) are important because they include the situation where / _1 (0) 
makes a non-crossing tangential intersection with the boundary of Bq. Detecting 
this situation is expensive, and provides a main motivation for exploiting isotopy. 
There is a similar expansion B + of B into B' for type (c), as indicated by dashed 
curves in Figure^c). In this case, we do not need to augment the graph G = (V, E) 
with any new vertices or edges, because G already contains an edge representing 
this component in B. 

We capture the above intuitive explanations as a lemma: 

Lemma 6. Assume that f is normalized relative to a quadtree T. Let B C Bq be a 
boundary box and let B' be its partner (a complementary box) that satisfies C\ or 
Cq. Then there exists a region B + such that: 

(i) B + is an expansion of B into B' : 

B C B+ C B U B' 

and the boundary of B + consists of two connected parts di(B + ) and d2(B + ): 
d\(B + ) is the union of three non-boundary sides of B and d2(B + ) is either the 
fourth side of B or a curve in the interior of B' (in Figure^ it is a dashed curve). 
In particular B+ n d(B') C B n B' . 

(ii) The augmented straightline graph G + = (V + ,E + ) when restricted to B is iso- 
topic to / _1 (0) n B + , i.e., 

G+ns«/- 1 (o)ns+. 

Moreover, the isotopy I on B + that witnesses this graph isotopy can be chosen to 
be the identity on di(B + ). 
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Proof. Note first that even though / has been normalized, it came from the 
situation where all boxes of T satisfy Co or C\. This, in turn, implies that the 
properties of Lemma [3] apply to incursions and excursions and that the isotopies 
chosen in the normalization process can be chosen to be the identity on these 
incursions and excursions. Note that G + f~l B has between and 3 components. 
This is because at most two edges can appear in GDB according to the Plantinga & 
Vegter construction rules, and in Types (a) and (b), the augmented graph G + has 
an additional component. We now consider each type in turn, and, in each case, 
define the expansion B + of B. These expansions are illustrated in Figure ^ a)- (e). 

For Type (e), only an excursion is possible. If there is an excursion, we can 
expand B into B + to ensure that / _1 (0) never intersects d 2 (B + ). 

For Types (c) and (d), the / _1 (0) intersects BOB' at least once. If it intersects 
BOB' multiple times, we can expand B into B + to ensure that / _1 (0) intersects 
the boundary d 2 (B + ) exactly once. 

For Types (a) and (b), there are two sides of B' where the curve / _1 (0) intersects 
at least once. Moreover, there is a unique connected component X of / _1 (0) n(B'U 
B) that connects these two sides (see Figure |6ja,b)). Recall that we say X has two 
possible dispositions: either X intersects the side B n B' or it does not. In either 
case, we can expand B into B + so that X intersects d2{B + ) exactly twice. This 
component X n B + is represented by the augmented edges (u,v), (v,w) in E + . 

Q.E.D. 

We now present the Extended Plantinga & Vegter algorithm. It has 3 

Phases that parallel the algorithm in Section 2. Phase i (for i = 1, 2,3) works off 
queues Qi and Q\, transferring boxes into Qi+\ and Q' i+1 . The queue Q\ holds 
complementary boxes while Qi holds regular boxes. Initially, Qi contains Bq and 
Q\ contains the four complementary boxes to Bq. 

• PHASE 1: SUBDIVISION. While Q x is non-empty, remove some B from 
Qi, and perform the following: If Cq(B) holds B is discarded. If C\{B) 
holds, and also Ci(B') or Co(B') holds for every complementary box B' 
of B, then (a) insert B into Q 2 , and (b) for each complementary box B' 
that satisfies C\ but not Co, insert B' into Q' 2 . Otherwise, split B into four 
subboxes which are inserted into Q\ and half-split B' (this means that we 
split it into four children and consider the two children which intersect 8Bq 
and put the two children into Qi). 

• PHASE 2: BALANCING. The balancing of boxes in Q 2 is done as in Phase 
2 of Section [3] Note that boxes are inserted into Q3 by this process. Next, 
we perform an analogous while-loop on Q' 2 : While Q' 2 is non-empty, remove 
any B' from Q' 2 . If B' and its partner are different sizes, we then half-split 
B' and put the children of B' into Q' 2 provided Co does not hold. Otherwise, 
we place B' into Q' 3 . 

• PHASE 3: CONSTRUCTION. First, perform Phase 3 of Section [3] which 
constructs a graph G = (V,E) from the boxes in queue Q3. Next we 
augment this graph using queue Q 3 : for each B' £ Q 3 , if B' is a transient 
type, we insert three vertices and two edges to the graph G as described 
above. The resulting straightline graph is denoted G + = (V + ,E + ). 
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4.3. Weak Correctness Statement. We are ready to prove the correctness of 
this Extended Plantinga & Vegter Algorithm. We define B' to be the union of B 
with all of the complementary boxes B' that were placed into Q' 3 . 

Before we formulate the correctness statement, we will discuss the implications 
of using an infinitesimal perturbation. In most cases the perturbation does not 
change the topology, e.g., when the curve intersects a corner in the interior of Bq, 
the topology in Bp do es not change. This is the reason that in Plantinga |2006 



Plantinga and Vegter 2004 , the curve S = / 1 (0) is isotopic to S = / 1 (0). In 
our setting, 5flBo may n0 longer be isotopic to S fl B . For instance, if S either 
makes a non-crossing tangential intersection with the boundary 8Bq at the corner 
of a box B of T or passes through a corner of Bq, this intersection is an isolated 
component of S(l Bq; we could lose or enlarge this component in SCiBq, depending 
on the sign of / nearby. To correctly analyze this case would require a more delicate 
consideration of the bounary, in particular, how complementary boxes of neighbors 
interact. Our correctness statement is therefore about SC\B , and not about SHBq. 

We regard the use of / as a reasonable compromise because (1) / is an infinites- 
imal perturbation of /; (2) / is easy to implement and is an effective perturbation 
(comparing favorably to the common alternative of, say, a randomized perturba- 
tion); (3) we know exactly when S fl Bq deviates from S H Bo (i.e., we encounter 
a zero at a box corner in the boundary of -Bo); (4) it is a very simple solution to 
what would otherwise be serious complications arising from various degeneracies 
(e.g., when S contains a horizontal or vertical component); and finally, (5) its use 
is consistent with the exploitation of isotopy in Plantinga & Vegter. 

Theorem 7 (Weak Correctness). Let the curve S = / _1 (0) be non-singular in the 
box Bq, and G + = (V + ,E + ) be the augmented straightline graph constructed by the 
Extended Plantinga & Vegter Algorithm. Then there exists a region Bq isotopic to 
Bo, 

Bo Q B^ C B'o, 

such that 

G + «r 1 (o)ns+. 

Proof. Let T be the quadtree at the end of the Extended Plantinga & Vegter 
algorithm and / be the standard perturbation of / relative to T. Using Theorem [i] 
we have 

(3) /- 1 (o)n J B «/r 1 (o)nBo, 

where /i is a normalization of / relative to T. Let Bq be the union of Bq with all of 
the caps constructed using Lemma [6j Since the isotopy constructed in Theorem [5] 
acts as the identity outside of Bo, it follows from ^ that: 

(4) r 1 (0)ni? o +«/r 1 (0)ni? o + . 

By a direct consequence of Theorem [5J we have G + (IB m / 1 ~ 1 (0) n B for each 
interior box B of T. For each boundary box B of T, Lemma [6] shows that G + flBsi 
/ 1 _1 (0) fl B + , where B + is the union of B and its cap. 

Now the composition of all the isotopies involved gives the isotopy G + ~ /-f 1 (0)n 
Bq, which, when combined with (W|), yields the desired result. Q.E.D. 
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REMARKS: 

1. Under the conditions of this theorem, we call G + a weak isotopic approxi- 
mation to / _1 (0) n Bq. In fact, once there is a collar around Bq where in each 
complementary box either Co or C\ holds, there is a spectrum of similar algorithms 
that allow for various correctness statements. This extension algorithm lies at one 
end of this spectrum, while the naive extension, which requires / (0) to intersect 
each boundary segment at most onces lies at the other end. The algorithms along 
this spectrum reflect a trade off between using isotopy on the boundary and topo- 
logical correctness within Bq. A slightly different approach, which is based on an 
algorithm that lies between these two extremes, was used in our ISSAC 2008 pro- 
ceedings version, that solution is in the middle of this spectrum: there, we required 
/ _1 (0) to only intersect OBq transversally. Then we can explicitly distinguish be- 
tween the two dispositions in the transient boxes (Types (a) or (b)). To do this, we 
require an iteration at each transient box. This application is a main reason why 
these boxes are called transient, after a finite number of additional subdivisions, the 
boxes will become non-transient. Finally the algorithm in ISSAC 2008 augments 
G with an edge if and only if this iteration detects an incursion. We prefer the 
current solution because it is non-trivial to ensure that / _1 (0) only intersects 8Bq 
transversally and because it uses isotopy as much as possible. 

2. We can make the collar B' Q \ B around B as narrow as desired: to ensure a 
perturbation bound of e > 0, it is sufficient that each of the complementary boxes 
where Co does not hold has width at most e/A. This is easily done by modifying 
the balancing phase of the algorithm. 

3. We have assumed that complementary boxes have the same width as their 
partner. An alternative, possibly more efficient, approach is to allow complemen- 
tary boxes to have widths less than their partners. Call these "subcomplementary 
boxes" . A boundary box B can have many subcomplementary boxes. We can do 
half-splits of subcomplementary boxes as long as their widths are greater than e. 
However, to ensure topological correctness, it becomes necessary to insist that a 
box has only a limited number of certain types of "subcomplementary boxes" . For 
instance, a boundary box B may be restricted to have at most one subcomplemen- 
tary box of Types (c) or (d) , and this occurs under strict conditions (we leave the 
details to the reader). 

4.4. Extension to Non-simply Connected Regions. It is essential in our ap- 
plications later to extend the above refinements to non-simply connected regions. 
Recall that we only consider regions i?o which come from subdivisions. To extend 
Theorem [7] to these regions, we note two simple modifications: 

(I) A complementary box B' of a boundary box B C R may intersect the inte- 
rior of Rq or other complementary boxes. Thus, Phase I must split such boundary 
boxes B sufficiently. Such interference checks can be checked during the subdivision 
phase. 

(II) The region Rq can have concave corners. A complementary box B' at a concave 
corner has two partners B C R Q . Relative to each partner B, we classify B' into one 
of 5 types as in Figure|6ja)-(e). This is illustrated in Figure [8fi)-(iv). In Figure[8]ji), 
for instance, the box B' is Type (a) (hence transient) relative to the indicated box 
B, but it is Type (d) (hence non-transient) relative to the other partner B. In Fig- 
ure |8](ii), the corner complementary box B' is Type (b) (hence transient) relative 
to both B and B. Similarly, the other two cases seen in Figure ISTiii)-(iv) have dual 
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classifications. We modify our augmentation of the graph G — (V, E) as follows: for 
each complementary box at a concave corner, we consider its classification relative 
to each choice of partner B: if the classification is transient, as before, we add three 
vertices u, v, w and edges (u, v), (v, w) to G on the side of B' n B. 



B 


B' 








'b ~ 


— 1 






(a)/(d) transient/non-trans. (b)/{b) transient (b)/(c) transient /non-trans. (c)/(c) non-trans 
W (ii) (ni) (vi) 



B 


H 7 — ' 

B' 


~* 





(e)/(e) non-trans. 
W 



Figure 8. Classification of complementary boxes at concave corners. 



5. Evaluation Bounds 
For any function /, define its evaluation bound to be 

EV(/):=inf{|/(p)| : f(p) ? 0, V/(p) = 0} 



Such bounds were used in [Cheng et al. 2007[ Burr et al. 2009a| . From Propo- 



sition [2j we see that {f(p) : V/(p) = 0} is a finite set and therefore EV(/) > 0. 
However there is no explicit bound readily available. The main objective of this 
section is to provide such a bound: 

Theorem 8. IffeZ[X,Y] has degree d and \\f\\ < 2 L then - lgEV(/) = O (<P(L+ 
dlogd)). More precisely, 

EV(Z)- 1 < max{[d 6 2 i + 2rf + 11 ]' i2 " 1 , {d 3d+s 2 3L+5d } d } 
Before giving the proof, we provide some definitions and preparation. Let h = 
V'/. , i: <i,,.\')' ' G C[X,Y}. \\h\\ k := yjEt+ J= oK\ k wiU denote the fc-norm of 
h, where we use k — 1,2. We just write \\h\\ for ||/i||oo := maxla^l, denoting the 
height of h. 

Now consider / as in the statement of the theorem considered as a function on 
C 2 . As input parameters in our bounds, we use d and L where deg/ < d and 
ll/H < 2 L . Let f x ,fy denote the derivatives of /. We may write 



hi 

ZERO(/j;, fy) = Ui U Vj 

i J 

where Ui are the 1-dimensional irreducible components and Vj are the 0-dimensional 
irreducible components. On each component Ui, one can show that the function / 
is constant. E.g. / = (xy + l) 2 — 1, f x = 2{xy + l)y and f y — 2{xy + l)x. Then 
U\ = {xy + 1 = 0} and V\ = {(0, 0)}. The function / is equal to —1 on U\ and 
on V\. 

Let g :— GCD(f x , f y ) and also 

fx-=fx/g, fv-=fv/g- 

Clearly, we have 

Zero(/ 2; , f y ) = Zero(5) U Zero(/z,/Q. 
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Since GCD(f x , f y ) — 1, we conclude that Zero(/ x , f y ) has no 1-dimensional compo- 
nents. Conversely, the hyper-surface ZERO(g) has no O-dimensional components, 
as a subvariety of C 2 . This proves: 

Lemma 9. 

Zero( 3 ) =\JUi, Zero(^, fy) =\JVj. 

i 3 

We now view the ring Z[X, Y] ~ Z[Y][X) ~ Z[A][F] in three alternative ways: 
A bivariate polynomial / in X and Y can be written as / = f(X, Y), f = f(Y: X) 
or / = f(X: Y) to indicate these three views respectively. As a member of Z[X, Y], 
the coefficients of /(X, Y) are elements of Z. But / = f(Y; X) is a member of 
Z[Y~][A] whose coefficients are elements of Z[Y}. The leading coefficient and degree 
of / are likewise affected by these views: lc(f{Y;X)) G Z[Y] but lc(f(X,Y)) G Z, 
e? = deg(/(A, F)) is the total degree of / while deg(f(Y;X)) is the largest power 
of X occurring in /. 

We use Mahler's basic inequality ( [Yap[ |2000] p. 351]) that if p G Z[A, Y] and 
p\q then 



(5) \\p(X,Y)\\ 1 <2 D \\ q (X,Y)\ 



where D = deg(q(X;Y)) + deg(q(Y;X)). This implies: 

(6) 11.9(^,^)11! < 4 rf - 1 d 3 2 i , 11^(^,^)11! < 4 d " 1 d 3 2 i . 

since 5 |/ x , and < d 2 ||/,|| < d 2 -d2 L , deg(f x (X; Y)) + deg(/ K (Y; X)) < 

2d — 2. The bound then follows from 

Let h{X) be the leading coefficient of g(X;Y). Since /i(A) has degree < d — 1, 
there is an integer a?o G {0, 1, . . . , d— 1} such that /i(xo) ^ 0. Intersect Zero(<?) with 
the line X = xq. We claim that this line cuts each non- vertical component Ui in a 
finite but non-zero number of points. In proof, let g = Y\ i g 3 i where ZERO(gi) = Ui. 
Setting di := deg gi{Y; X), we see that the vertical components correspond to = 0. 
Then lc(g(X;Y)) = l\ i lc{g l (X;Y)) and h(x ) = lc(g(x ;Y)) ^ iff for all i, 
lc(gi(xo; Y)) ^ 0. But gi(xo; Y) is a polynomial of degree di in Z[Y], and, therefore, 
has exactly di solutions in C. 

Write f (Y) :=f(x ,Y) and g (Y) = g{x .Y). From (§: 

lbo||i <d d \\g(X,Y)\\i <4 d - 1 d d+3 2 L . 
It is also easy to see that 

ll/oll < d d+1 2 L+1 . 

Suppose f3 G Zero (50) \ Zero(/ ). We want a lower bound on |/q(/3)|. For this 



purpose, we use an evaluation bound from [Burr et al. 2009a Theorem 13(b) 



Proposition 10 (Evaluation Bound Burr et al. 2009a ). Let <f>{x),T]{x) G C[x] be 



complex polynomials of degrees m and n. Let f3\,. . . ,f3 n be all the zeros of r)(x). 
Suppose there exists relatively prime F, H G Z[x] such that F = cjxp, H = rffj for 
some cj>, rj G C [x] . If the degrees of <j> o,nd rj are rh and n, then 
n 1 
E[l^)l> 



~ / „ \ n+n 

lc(rj) m ■ ((m + l)||0||)™ M(rj) m ■ Km + l)||^||J M{H) 
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Here the Mahler measure M(h) of a function h 6 C[x] with zeros a%, ■ ■ ■ , a n , is 
defined as M(h) :— \lc(h)\ II|a |>i \ a i\- We shall choose the functions in Proposi- 
tion [TO] as follows: 

H:= GCD(/o, 9o )- 
Moreover, let <fi := 1, := Y — f3 and 77 := Hjr\ 6 C[x]. Hence 
m < d, n = 1, m = 0, n < d — 1. 

Also 

lc{rj) = lc(H)=lc(g ) < \\g \\ < \\g \\i- 

Further, 

M(5j)<A/(ff)<||fl-||i< 2 d -|| 5o ||i. 
Finally, an application of Proposition [lO] gives 

l/oC^r 1 < lc(7i) d -((d+l)\\f \\) d - 1 -M(ri) d 

< Mv) ■ d\\M\ ■ M(?j)} d ( as (d + < d d for d > 2) 

< [!| 50 ||i-drf d+1 2 L - 2 d || ff o||i] d 
(7) < [d 3d + 8 2 3L+5d ] d . 

Q is a lower bound on |/(p)| where p lies in a non- vertical component Ui. By 
considering g(Y;X), the same bound applies for |/(p)| when p lies in a vertical 
component Ui. 

Now, we obtain a lower bound for f(p) with p G Zero(/ x , Consider the 
system E C Z[X, Y, Z] where 

s = {z - /(x, y), 7^(x, r)} 

The zeros (6,6,6) = (x,y,f(x,y)) £ C 3 of E satisfy 6 = /(6,6)- Since S is a 



zero dimensional system, we may apply the multivariate zero bound in Yap 2000 
p. 350]. This bound says that 

< {2 3 ' 2 NK) D 2^ d -^ 
where N = ( 1+2( 3 d_1) ), D = d 2 - 1 and 

K = max{V3, ||^|| 2 , ||^|| 2 , \\Z - f(X, Y)\\ 2 }. 

We have \\Z - f{X.Y)\\ 2 < (d+l)2 L+1 . From |6j), we see that K < A d d 3 2 L . Using 
the bound N < 2d 3 , we obtain 

(8) < [2 2 • 2d 3 ■ 4 d d 3 2 L f^ ■ 2^-^ < [d^+^+nf-i^ 

Now Theorem [8] easily follows from Q and pL 

6. Isolating Singular Points 

In the remainder of this paper, we assume that / € Y] and allow the curve 
S = / _1 (0) to have singular points, except on the boundary dB . We would like 
to use the Extended Plantinga & Vegter algorithm to compute an isotopic approx- 
imation to Zero(/) when / has only isolated singularities. Since the Plantinga & 
Vegter algorithm does not terminate near singular points, it is necessary to isolate 
the singular points from the rest of Bq. 

We use the auxiliary function F = f 2 + f% + f 2 - Finding the singular points of 
/ _1 (0) amounts to locating and isolating the zeros of this non- negative function. 
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We use a simple mountain pass theorem Jabri 2003 adapted to Bq to ensure our 
algorithm isolates the zeros. 

Theorem 11. Suppose that F > on Bq, and that F > on dB a . Then for any 
two distinct roots p,q of F in Bq, there exists a continuous path 7 : [0,1] — > Bq 
connecting p and q which satisfies the following: 

• 7 minimizes M 7 := max xg jo,i] F(l( x )) among all paths connecting p and q 
in Bq. 

• 7 contains a point y such that either VF(-y(y)) = or 7(2/) € dB. 

This can be proved using path deformation and the compactness of Bq, or it can 
be seen as a simple application of the topological mountain pass theorem presented 
in Jabri 2003 . Because of this theorem, distinct zeros of F within Bq are separated 
by barriers of height e = min(EV(F), minF(9i?o))- This leads us to the following 
multistep process to localize these zeros. The goal is to hnd a small rectangle with 
diameter less than some S around each zero. 

STEP 0: DETERMINING e. Initialize e to any lower bound on EV(F). Also, 
initialize Qq to be {Bq} and Q\ to be empty. Qq will contain boxes S that intersect 
dB and have £ DF(S); these boxes will then be subdivided in order to compute 
a lower bound on F(dB ). While Qq is non-empty, remove a square S from it 
and evaluate □-F(S'). If 0F(S) > we push S into queue Q\ and also update e to 
min{e, minDF(S')}. If € QF(S), subdivide S and push the children of S which 
intersect 8Bq into Qq, and the others into Q\. When Qq is empty, we stop and fix 
the current value of e for the remainder of the algorithm. 

STEP 1: INITIAL SUBDIVISION. Initialize queue Q 2 to be empty. While there 
is an S in Q 1 , remove it and evaluate QF(S). If UF(S) > e/2, then discard S. Else 
if □-F(S') < e, push S into Q 2 . Else subdivide S and push its children into Q%. 

Once Q\ is empty, group the elements of Q 2 into connected regions A4 (i € 
7). Each Ai contains at most one root, since otherwise, there would be a path 
connecting the roots within A- The value of F along this path would be less than 
e, contradicting the mountain pass theorem. Let C be the region B \ UjAj- Note 
by Step that F is greater than e/2 on C and that 8Bq C C, i.e., each Ai doesn't 
intersect 8B . 

STEP 2: REFINEMENT. For each A, (i e I), initialize queue Q 2 a with all 
squares S € Aj. So long as neither terminating condition 1 nor 2 (below) hold, 
we perform the following: For each S in Q 2 ,i, if G DF(S), subdivide S and push 
its children into Q 2 ,i- If ^ QF(S), discard S. We terminate when either of the 
following two conditions are met: 

(1) Qi.i is empty, in which case there isn't a zero in A4. 

(2) A\, the contents of Q 2 ,i satisfy all of the following: 

(a) UF(S) < e/2 for some S € A\ 

(b) Ri, the smallest rectangle containing A\, lies within the region covered 
by the original Ai. 

(c) The diameter of Ri is less than 8. 

It is clear from the definition of F that this step will halt. We claim that each 
Ri contains exactly one root. In Step 1, we showed that A, contains at most one 
root. To see that Ri contains a root, take a point of A\ where e/2, then follow the 
path of steepest descent to reach a zero of F. Because F is less than e/2 on this 
curve, the curve cannot pass through the region C to reach any other Rj or to leave 
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Bq. Therefore there must be a zero within Ai. It is in Ri because our conditions 
ensure that F is positive on A t \ Ri. 



7. Determining the Degree of Singular Points 



The following standard result from Krantz and Parks 1992 Lojasiewicz 1991 



shows that the global structure of zero sets: 

PROPOSITION 12 (Zero Structure). Let f be real analytic. Then Zero(/) can be de- 
composed into a finite union of pieces homeomorphic to (0, 1), pieces homeomorphic 
to S 1 , and singular points. 

In our current situation, the pieces which are homeomorphic to (0, 1) are smooth 
open subsets of the irreducible components of Zero(/). 

Viewing Zero(/) as a multigraph H, the branching degree of a singular point is 
its degree as a vertex of H. We now determine such degrees. Let 83 be a separation 
bound between singular points, so if p and q are two distinct singular points of 
Zero(/), then the distance between p and q is at least S3. Let ^4 be a separation 
bound so that if r is a point on Zero(/) such that V/(r) is parallel to the line 
between r and a singular point p, then the distance between p and r is at least 
54. If s is on Zero(/) so that the distance between s and a singular point p is 
smaller than either £3 or 84, then by following the paths Zero(/) away from s, one 
of the paths strictly monotonically approaches p until it reaches p and the other 
path locally strictly monotonically recedes from p. |Yap , 2006] provided an explicit 



bound on S3 as a function of degree and height of f(X, Y): 

83 > min{(16 d+2 256 i 81 2d d 5 )- d ,(2 8i+21 3 8d )- 2 } . 

7.1. Lower Bound on ^4. To derive an explicit bound on £4, we consider the 
following 6 polynomials in Z[z,p x ,p y , q x , q y \: 

(9) I f :={/(p), /(?), f x (q), f y (q), (p - q) x f v (p) - (p - q)yf x (p),Z 2 - \\p - qf} 

where p = (p x ,p y ) and q = (q x ,q y ) are points on Zero(/). q is a singular point 
of Zero(J). Moreover, define "(p — q) x " to mean p x — q x , so that the equation 
(P - q)xf y (p) ~ (p ~ q) v fx{p) = implies that (p - q) is parallel to V/(p). The 
equation z 2 — \\p — q\\ = implies that z is the distance between p and q. 

Consider the projection IL[Zero(//)] of the zeros of // onto the z-coordinate. 
Then is inf{|z| : z £ IT z [Zero(//)], z ^ 0}. We obtain a lower bound on S4 using 
the following general theorem from |Brownawell and Yap 2009 : 



PROPOSITION 13. Let I:=(P 1} , . . . , , P m ) C A:=Z[Xx, . . . ,X n ]. Let *P be an iso- 
lated prime component of I whose projection onto the first coordinate, IIi(Zero(^P)), 
is a finite set. If £ = (d, ,. . . , , £„) € Zero(*}3) and £i ^ 0, then 

I Cx I > ((n + l)2 e n+2)-"(«+l)D"- fc 

where 

• dim*P = k. 

• H > Height (Pi), and 

• D > deg(Pi), i=l, m. 
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If II z [Zero(//)] is a finite set, then we use the bounds n — 5, m — 6, k < 4, 
D = max{2, d}, and H < d2 L+1 in applying this theorem to get the following 
bound: 

S 4 > (6 2 e 7 )- 30Z55 (4 4 - 6-d2 i+1 )- 5D4 
Combining the two cases for the value of D gives 

<5 4 > min {(6 2 e 7 )- 30dB (4 4 • 6 • d2 L+1 )~ 5d \ (6 2 e 7 r 30 ' 25 (4 4 ■ 6 • d2 i+1 )- 5 ' 24 } 

It remains to show that II 2 [Zero(//)] is a finite set. We prove this in the following 
lemma: 

Lemma 14. IT z [Zero(//)] is a finite set. 

Proof. Without loss of generality, we apply a translation so that we can assume 
that q is at the origin. To show that this image is a finite set, we show that 
Zero(/ (p) , p x f y (p) —Pyfx(p)) is contained in finitely many circles centered at the 
origin. Then, the possible values of z are the radii of these circles, of which there 
are finitely many. 



By Proposition 12 we know that each component of Zero(/) is either a single 
point or a smooth one dimensional manifold without boundary. Since there are 
finitely many components which are a single point, these components are contained 
within finitely many circles centered at the origin. 

Take any one dimensional component M and let r, s be two points in M and 
7 : [0, 1] — > M a smooth path from 7(0) = r to 7(1) = s in M . Also, define 
p : C 2 —> C by p(x,y) = x 2 + y 2 . For any p £ M, we note that since M is 
smooth, the tangent line of M (or equivalently Zero(/)) at p is perpendicular to 
V(/) = {fx(p), fy{p)) and this gradient is non-zero because p is a smooth point of 
Zero(/). In addition, since Pxf y (p)-p y fx(p), we know that (jp x ,p v ) || (fx(p), /(?/))• 
Also note that Vp = (2x, 2y) = 2 ■ (x, y). 

Now, we consider the square of the distance between 7 and the origin, pip/it)). 
Taking the derivative of this function gives 

j t Ph(t)) - (Vp)( 7 (t)) ■ V(t) = 2 7 « • i{t) = 0. 

Since this derivative is always zero, it implies the distance from the origin to points 
7(t) on M is constant. Therefore, M is contained in a unique circle centered at the 
origin. Q.E.D. 

To find the degree of a singular point, assume that we have two boxes E>i -2 B 2 
where the diameter of B\ is less than both ^3 and B2 contains a singular point 
of / and there is some radius r > such that a circle of radius r centered at any 
point p inside B2 must lie entirely within the annulus B\ \ B2, see Figure |9ja). 
Note this condition is satisfied if B\ is at least 5 times larger than B2 and B2 is 
centered in Bi, which is the typical situation we consider below. See Figure J9jb). 
Furthermore, to apply our extended Plantinga & Vegter algorithm of Section 4, we 
ensure that B\ \ B2 comes from a subdivision. 

Now, there are 3 types of components (other than isolated points) in Zero(/) n 
(B\ \ int(i?2)): (1) images of [0, 1] both of whose endpoints are on dB\, (2) images 
of [0,1] both of whose endpoints are on dB 2 , and (3) images of [0,1] with one 
endpoint on each of dBi and dB 2 - These three types are illustrated in Figure |9ja). 
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Figure 9. Annular region B\ \ Bi with singularity p and the three 
types (1), (2), (3) of components. 



Let s be a point on any of these components, then traveling along Zero(/) in 
one direction must lead to the singular point and the other direction must leave the 
neighborhood (be further than min{<!>3, £4}) of the singular point. For, if not, then 
there must be a point r on Zero(/) (1 B\ such that V/(r) is in the same direction 
as the line between s and the singular point, which is impossible since the width 
of B\ is smaller than 84. Now, any piece defined by Proposition [T2| which reaches 
the singular point exits the neighborhood of the singular point and the only way 
to leave the neighborhood is by way of a type (3) component. In addition, each 
piece includes either one or two components of type (3), and it only includes two 
components if both endpoints reach the singular point. This follows from the choice 
of 84. This shows: 

Lemma 15. The degree of the singular point in B 2 is the number of components of 
type 3. 

8. Overall Algorithm 

We now put all the above elements together to find a weak isotopic approx- 
imation to the algebraic curve S = f~ l (0) within a region R , coming from a 
subdivision, where f(X,Y) € Z[X,Y] has only isolated singularities. We first find 
the singularities of the curve S in i?o- Using the technique of Section 5, we can 
isolate the singularities pi (i — 1,2,...) into disjoint boxes Bi. We assume the 
width of the Bi's is at most min{($3, 84} /6. Let B[ be the box of width 5 times the 
width of Bi, and concentric with Bi\ we further assume B[ C R . Note that these 
combinations are chosen to ensure that we have the typical situation in Section [7] 
Now we proceed to run the extended Plantinga & Vegter algorithm on the nice 
region R* := R \ \J t Bi, yielding a polygonal approximation G. We directly incor- 
porate the technique of Section [7] into the following argument. If pi is the singular 
point in Bi, then the degree of pi is equal to the number of type (3) components in 
G n {B[ \ Bi). We connect these components directly to p i: and discard any type 
(2) components. This produces the desired isotopic approximation. 

Remarks: We have not discussed e-approximation because this is relatively easy 
to achieve in the Plantinga & Vegter approach. We only have to make sure that 
each subdivision box that contains a portion of the polygonal approximation G 
has width at most e/4 since the result of the Plantinga & Vegter algorithm only 
deforms the original curve at most one cell away. 
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9. Conclusion 

This paper presents the first complete numerical subdivision algorithm for mesh- 
ing an implicit algebraic curve that has only isolated singularities. This solves an 



open problem in the exact numerical approaches to meshing in 2-D Boissonnat 



et al. 2006 p. 187]. We pose three challenges: 
(a) A worst case complexity bound for our procedure is possible. But this may not 
be the best way to measure adaptive algorithms. We would like to provide adaptive 



bounds, similar to the integral analysis in |Burr et al.| 2009a| for 1-D problems. 



(b) In 3-D, a square-free integer polynomial f(X,Y,Z) could have a 1-dimensional 
singularities. We pose the problem of designing a purely numerical subdivision al- 
gorithm to handle 1-dimensional singularities. 

(c) The practical implementation of an adaptive algorithm handling singularities, 
even based on our outline, must handle many important details. Computational 
experience is invaluable for future research into singularity computation. 



References 

Saugata Basu, Richard Pollack, and Marie-Frangoise Roy. Algorithms in Real Alge- 
braic Geometry. Algorithms and Computation in Mathematics. Springer, 2003. 

J.-D. Boissonnat and S. Oudot. Provably good sampling and meshing of Lipschitz 
surfaces. In Proc. 22nd ACM Symp. on Comp. Geometry, pages 337-346, 2006. 
Sedona, Arizona. 

J.-D. Boissonnat, D. Cohen-Steiner, B. Mourrain, G. Rote, and G. Vegter. Meshing 
of surfaces. In J.-D. Boissonnat and M. Teillaud, editors, Effective Computational 
Geometry for Curves and Surfaces. Springer, 2006. Chapter 5. 

Jean-Daniel Boissonnat, David Cohen-Steiner, and Gert Vegter. Isotopic implicit 
surfaces meshing. In ACM Symp. Theory of Comput., pages 301-309, 2004. 

W. D. Brownawell and Chee K. Yap. Lower bounds for zero-dimensional projec- 
tions. In Proc. 34th Int'l Symp. Symbolic and Algebraic Comp. (ISSAC'09), pages 
79-86, 2009. KIAS, Seoul, Korea, Jul 28-31, 2009. Submitted, JSC. 

Michael Burr, Felix Krahmer, and Chee Yap. Continuous amortization: A non- 
probabilistic adaptive analysis technique. Electronic Colloquium on Computa- 
tional Complexity (ECCC), TR09(136), December 2009a. URL |http://eccc.| 
|hpi-web . de/report/2009/136/[ 

Michael Burr, Vikram Sharma, and Chee Yap. Evaluation-based root isolation, 
February 2009b. In preparation. 

Jin-San Cheng, Xiao-Shan Gao, and Chee K. Yap. Complete numerical isolation 
of real zeros in general triangular systems. In Proc. Int'l Symp. Symbolic and 
Algebraic Comp. (ISSAC'07), pages 92-99, 2007. Waterloo, Canada, Jul 29- Aug 
1, 2007. DOI: http://doi.acm.org/10.1145/1277548.1277562. In press, Journal of 
Symbolic Computation. 

S.-W. Cheng, T.K. Dey, E.A. Ramos, and T. Ray. Sampling and meshing a surface 
with guaranteed topology and geometry. In Proc. 20th ACM Symp. on Comp. 
Geometry, pages 280-289, 2004. 

D. Cox, J. Little, and D. O'Shea. Ideals, Varieties and Algorithms: An Introduc- 
tion to Computational Algebraic Geometry and Commutative Algebra. Springer- 
Verlag, New York, 1992. 



26 



MICHAEL BURR, SUNG WOO CHOI, BEN GALEHOUSE, AND CHEE YAP 



Zilin Du and Chee Yap. Uniform complexity of approximating hypergeometric 

functions with absolute error. In Sung-il Pae and Hyungju Park, editors, Proc. 7th 

Asian Symp. on Computer Math. (ASCM 2005), pages 246-249, 2006. 
Joe Harris. Algebraic Geometry. Springer- Verlag, New York, 1992. 
Robin Hartshorne. Algebraic Geometry. Springer- Verlag, New York, 1977. 
H. Hong. An efficient method for analyzing the topology of plane real algebraic 

curves. Mathematics and Computers in Simulation, 42:571-582, 1996. 
Youssef Jabri. The Mountain Pass Theorem: Variants, Generalizations and Some 

Applications. Encyclopedia of Mathematics and its Applications. Cambridge 

University Press, 2003. 
S. G. Krantz and H. R. Parks. A primer of real analytic functions. Birkhauscr 

Verlag, Basel, 1992. ISBN 3-7643-2768-5. 
Long Lin and Chee Yap. Adaptive isotopic approximation of nonsingular curves: 

the parametrizability and non-local isotopy approach. In Proc. 25th ACM Symp. 

on Comp. Geometry, pages 351-360, June 2009. Aarhus, Denmark, Jun 8-10, 

2009. Accepted for Special Issue of SoCG 2009 in DCG. 
S. Lojasiewicz. Introduction to complex analytic geometry. Birkhauser Verlag, Basel, 

1991. ISBN 3-7643-1935-6. Translated from the Polish by Maciej Klimek. 
W. E. Lorensen and H. E. Cline. Marching cubes: A high resolution 3D surface 

construction algorithm. In Maureen C. Stone, editor, Computer Graphics (SIG- 

GRAPH '87 Proceedings), volume 21, pages 163-169, July 1987. 
Bernard Mourrain and J. -P. Tecourt. Isotopic meshing of a real algebraic surface. 

Technical Report RR-5508, INRIA, Sophia-Antipolis, France, February 2005. 

Also, electronic proceedings, MEGA 2005. 
Simon Plantinga. Certified Algorithms for Implicit Surfaces. Ph.D. thesis, Gronin- 

gen University, Institute for Mathematics and Computing Science, Groningen , 

Netherlands, December 2006. 
Simon Plantinga and Gert Vegter. Isotopic approximation of implicit curves and 

surfaces. In Proc. Eurographics Symposium on Geometry Processing, pages 245- 

254, New York, 2004. ACM Press. 
Elmar Schoemer and Nicola Wolpcrt. An exact and efficient approach for computing 

a cell in an arrangement of quadrics. Comput. Geometry: Theory and Appl, 33: 

65-97, 2006. 

Raimund Seidel and Nicola Wolpert. On the exact computation of the topology 
of real algebraic curves. In Proc. 21st ACM Symp. on Comp. Geometry, pages 
107-116, 2005. Pisa, Italy. 

J. M. Snyder. Generative Modeling for Computer Graphics and CAD: Symbolic 
Shape Design using Interval Analysis. Academic Press, 1992a. 

J. M. Snyder. Interval analysis for computer graphics. SIGGRAPH Com- 
put.Graphics, 26(2):121-130, 1992b. 

Barton T. Stander and John C. Hart. Guaranteeing the topology of an implicit sur- 
face polygonalization for interactive meshing. In Proc. 24th Computer Graphics 
and Interactive Techniques, pages 279-286, 1997. 

Chee K. Yap. Fundamental Problems of Algorithmic Algebra. Oxford University 
Press, 2000. 

Chee K. Yap. Complete subdivision algorithms, I: Intersection of Bezier curves. In 
22nd ACM Symp. on Comp. Geometry, pages 217-226, July 2006. 



ISOTOPIC MESHING OF SINGULAR ALGEBRAIC CURVES 



27 



Courant Institute, NYU 251 Mercer Street, NY, NY 10012. 

Current address: Mathematics Department, Fordham University, Bronx, NY 10458. 

E-mail address: burr@cims.nyu.edu 

Department of Mathematics, Duksung Women's University, Seoul 132-714, Korea. 
E-mail address: swchoi@ducksung.ac.kr 

Courant Institute, NYU 251 Mercer Street, NY, NY 10012. 

Current address: Max-Planck-Institut fiir Informtik Department 1: Algorithms and Complex- 
ity, 66123 Saarbriicken, Germany. 

E-mail address: bgalehou@mpi-inf.mpg.de 



Korea Institute of Advanced Study, Seoul, Korea and Courant Institute, NYU 251 
Mercer Street, NY, NY 10012. 
E-mail address: yap@cs.nyu.edu 



